Connected network of minima as a model glass: long time dynamics 



L. Angelani 1 , G. Parisi 2 , G. Ruocco 1 and G. Viliani 3 
1 Universitd di L'Aquila and Istituto Nazionale di Fisica della Materia, 1-67100, L'Aquila, Italy. 
Universitd di Roma La Sapienza and Istituto Nazion* 
3 Universitd di Trento and Istituto Nazionale 



oo 

On 
<3\ 



(N 



I 

C/3 



I 

C 

o 
o 



> 

in 

m 
O 
oo 

On 
c3 



i 

a 

o 

o 



X 
S3 



e Fisica Nucleare, 1-00185, Roma, Italy. 
Fisica della Materia, 1-38050, Povo, Trento, Italy. 
(February 1, 2008) 



A simple model to investigate the long time dynamics of glass-formers is presented and applied to 
study a Lennard- Jones system in supercooled and glassy phases. According to our model, the point 
representing the system in the configurational phase space performs harmonic vibrations around 
(and activated jumps between) minima pertaining to a connected network. Exploiting the model, in 
agreement with the experimental results, we find evidence for: i) stretched relaxational dynamics; 
ii) a strong T-dependence of the stretching parameter; Hi) breakdown of the Stokes- Einstein law. 



PACS Numbers : 61.20.Lc, 64.70.Pf, 82.20.Wt 

In recent years many efforts were devoted to the un- 
derstanding of the phase space landscape in supercooled 
liquids and structural glasses, and, in particular, to the 
identification of those landscape details that are respon- 
sible for the structural arrest taking place at the glass 
transition temperature, T g GJ|. It has been recently spec- 
ulated pi that the free energy landscape of structural 
glasses is similar to that of some generalized spin glasses 
models, where it was shown that exists a dynamical tem- 
perature To (which is well defined in mean field approx- 
imation and becomes a crossover region in real systems) 
below which the dynamics is dominated by long time ac- 
tivated processes consisting of jumps among different free 
energy minima. The parallel between To and the criti- 
cal temperaure, Tc, of the mode cupling theory (MCT) 
H is natural. If verified, this parallel could lead to a 
microscopic description of the dynamic slowing down, 
predicted by MCT, in term of free energy landscape. 
Approaching Tc (or To), however, the presence of this 
very slow dynamics makes the numerical investigation of 
structural glasses very hard. To our knowledge, only few 
attempts were made in this direction pj. 

In this letter we introduce a new method to study 
the slow dynamics in glasses and in deeply supercooled 
liquids; at variance with the usual Molecular Dynamics 
(MD) simulation we describe the dynamics of the sys- 
tem as relaxations taking place in a connected network 
of potential energy minima. The jumps among minima 
are described by an appropiate master equation, and, in 
this way, we can investigate the long time behaviour of 
a glass in short simulation times as the solution of the 
master equation is an eigenvalue problem. The character- 
istics and connectivity of the minima, and other energy- 
landscape properties entering the determination of the 
transition probabilities, are inferred from the MD inves- 
tigation of a small system (one component Lennard- Jones 
in the present case). 

The physical quantities (total energy, pressure, trans- 
port coefficent) obtained from the model agree with those 
derived from MD up to temperature above the melt- 



ing point, supporting the jump model even in the liquid 
phase of Lennard- Jones fluids. In the low temperaure 
region we found evidence for: i) stretched behaviour of 
the relaxation process; ii) temperature dependence of the 
stretching exponent, (3k , which changes from w 1 at high 
T down to Ri 0.3 at low T; and Hi) breakdown of the 
Stokes-Einstein relation. All these results are in agree- 
ment with the experimental findings in real "fragile" 0] 
glasses. 

In a glass, the atoms are (almost) frozen in some 
(meta)stable positions. The short-time dynamics is dom- 
inated by small vibrations around the stable position. 
This dynamics can be described within the harmonic ap- 
proximation, and all the relevant information is obtained 
by diagonalizing the dynamical matrix. At long time, 
collective jumps among different stable positions involv- 
ing many atoms become possible and are controlled by 
a master equation. These long time relaxations are only 
apparently in contraddiction with the harmonic vibra- 
tional dynamics within the local minima; indeed, it has 
been recently shown that the relevant classical path for 
relaxation between adjacent minima is practically decou- 
pled from the other degrees of freedom, which are al- 
most harmonic ||. The transition rates, in turn, are de- 
termined by minima energies, barrier heights and other 
topological properties. 

In order to set up the connected netwok of minima and 
to determine the transition rates we need the topology 
of the multidimensional potential energy hypersurface of 
the system. To this end, we numerically analyse the Po- 
tential Energy Landscape (PEL) of small (N = 11 -j- 29 
atoms) Lennard- Jones systems whit periodic boundary 
conditions. The small size of the system allows us an ex- 
haustive investigation of the landscape but, at the same 
time, exhibits complex enough features and behaviour 
to capture the physics of the system. The atoms in- 
teract via the 6—12 Lennard- Jones potential Vlj{t) = 
4e[(f ) 12 - (s ) 6 ], whit e/K B = 125.2 K (K B Boltzmann's 
constant) and a = 0.3405 nm, appropriate for Argon. 
The simulated density is p — 42 mol/dm 3 . Due to the 
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small dimensions of the system, we choose the multi- 
image method, namely each particle interacts with many 
images of each other (in practice, due to the cutoff dis- 
tance that we impose on the pair potential, r cut = 2.6 a, 
there are at most 3 interacting images). 

The first step is to search for the potential energy min- 
ima. We perform a modified steepest descent method 
procedure , starting from high temperature MD config- 
urations to find the inherent configurations correspond- 
ing to local minima which often result to be crystalline- 
like. In order to establish whether a minimum corre- 
sponds to a glassy structure, we use the static stucture 
factor S(q) — N^ 1 | e'' Tl | 2 . For a pure crystalline 
configuration of N particles S(q) is made up by 'Bragg' 
peaks and its value at the peaks is S max = N, whereas 
for a glass one usually finds S max « 2 -j- 3. In small 
size samples there are obviously intermediate situations 
and for a minimum to be 'glassy' we adopt the criterion 

Smax < N/2. 

As a second step, for each pair of minima, a and b, 
with energy E a and E respectively, we first determine 
the mutual distance d a .b = min(\r a — r b \), where r a is the 
position vector of the minima in the 3iV dimension con- 
figurational space, and min indicates the minimization 
with respect to the all symmetry operations: continuous 
translations, permutations of particles and the 48 sym- 
metry operations of the cubic group. Then we analyse 
the potential energy profile experienced by the system 
in travelling from one minimum to another and deter- 
mine the potential energy barrier. Among the different 
paths joining a and b, we assume Q that the system 
follows that with the least action. The action integral 
is defined as S(£) = J, ds^V^s)) — V , where £ indi- 
cates a generic path, s the curvilinear coordinate, and 
V = min{i? a , Eb}. The minimization of S(£) is per- 
formed by dividing the path in n — 16 intervals and min- 
imizing the action function with respect to the extrema 
of the n segments constrained to move in iperplanes per- 
pendicular to the straigh path. The highest energy value, 
V a b, along the least action path (LAP) determines the 
saddle point of the path. Not all the pairs of minima 
are directly connected, since, sometimes, the LAP join- 
ing them crosses a third minimum. Therefore, there is 
a non-trivially connected network of minima. Next we 
mesure the curvature, defined as the determinant of the 
Hessian of potential energy function, in each minimum 
a (det{V a "}), and a - b saddle point (det{V;'£}). Also 
important is the absolute value of the negative curvature 
on the saddle point, Co a b- 

In order to give a full statistical description of the PEL, 
we study the distributions, P(x), of the relevant quan- 
tities x (here x represents E a , AE a b — E a — Eb, d a b, 
V a b, det{V"}, detjV^}, or ui a b) and their cross correla- 
tions, P(xi,X2)- Cross-correlations among the measured 
quantities are observed. The most evident is a linear 
correlation in double log scale between the distance and 



the barriers' height along the LAP between two minima. 
A rather weak correlation is also observed between the 
energy and curvature of extrema points of PEL. 

The model we introduce is a connected network of po- 
tential energy minima with a jump-dynamics described 
by an appropriate master equation: 



p a (t) = X c W ac p c (t) 



(1) 



where p a (t) is the probability that the system is in min- 
imum a at time t (actually p a (t) = p a (t\b), indicat- 
ing that at t = the system was in minimum b) and 
the non-diagonal elements of the transition matrix, the 
transition rates W a b, are determined from the energetic 
and topological properties of the PEL. In order to sat- 
isfy the equilibrium condition, p° — p a (t — > oo) oc 
(det V;")" 1/2 e-P E % with (3 = {KbT)' 1 , W ab must obey 
the detailed balance W a b Pi — Wb a P° a - Following [0 we 
make the ansatz: 



W ab = -2* 

7 



det V b " 
dct VK 



1/2 



,-P(E ab -E b ) 



(2) 



where 7 is a friction constant which actually determines 
the time scale. This choice of the transition matrix is 
based on the approximation of the problem of escape 
from a metastable state as a Markovian Brownian mul- 
tidimensional motion in the overdamped friction regime 

To set up our model minima-network we proceed in 
the following way. Having fixed the number of minima 
(M=400 in the present case), we extract the energy of 
these minima and their curvatures from the previously 
found bivariate distribution. For each minimum we ran- 
domly extract 20 minima connected to it and define a 
connection matrix c a b that contains the number of steps 
required to go from a to b. We then define the distance 
d a b as c a b times the value extracted from the distribu- 
tion of the distances P(d a b), and from these the ener- 
gies of saddle points. The further statistical features of 
saddle points (curvatures) and minima (trasverse com- 
ponent of the microscopic stress tensor, see below) are 
determined from bivariate (correlation curvature-saddle 
point energy) and simple extractions, respectively. 

To check the reliability of the model, we first concen- 
trate on the static properties. Following our model, the 
partition function is approximated by a sum over the 
minima and the harmonic vibrations around them: 



z{0) = /r 3iV/2 £ a (detv;")- 1/2 e"^ 



(3) 



In fig. I we show the potential energy of a LJ sys- 
tem with N — 29 particles as obtained through MD and 
as calculated from (||) by taking into account either all 
minima (dotted line), or only the glassy ones (full line). 
The MD data are obtained progressively heating the glass 
(o) up to the liquid phase, and then cooling the system 



2 



slowly (•), to obtain crystallization. We observe a quan- 
titative agreement between the MD data and the model 
up to T » 150 K, a temperature well above the melting 
point (T m « 80 K). At higher temperature the simple 
local- vibration/collective-jumps model fails. 

As for the dynamical properties, the solution of the 
master equation is found by numerical knowledge of 
eigenvalues and eigenvectors of the transition matrix X(n) 
and v a (n), with n = 1...M. In particular: 

Pa(t\b)=Y, Vain) o bin) Z X{n)t - (4) 

It is then possible to determine the statistical equilibrium 
average of a generic observable 0(t) from the knowledge 
of its value, O a b, calculated at the minima a and b: 

b a 

We consider three observables: the mass diffusion coef- 
ficient D, the shear viscosity n, and the structural relax- 
ation time r. In the first case O a b =\ L a ~Lb I ■> an d D = 
limt^oo < 0(t) > /6t. For the two other cases, we first 
determine the correlation function, of the off-diagonal el- 
ements of the stress tensor, C(t) =< <r zx (0)<t zx (t) >, 
with |]: 

^=-£^n^-)- (6) 

In the notation of eq. (||), O a b = fffaj 1 , (where a a is 
the value of the stress tensor calculated at minumum a) . 
Then the shear viscosity is calculated as 

/>oo 

n = (KbTV)- 1 / dt c(t), (7) 

Jo 

and the relaxation time r is derived from a fit of C(t) to 
a stretched exponential decay C(t) — C(0) exp — (t/r) l3K . 

We report the values obtained for M — 400 and 
avereged over 50 different extraction of the network pa- 
rameters. In fig. 2 we show the normalized correla- 
tion functions C(t)/C(0) calculated at different temper- 
atures together with their best fits. In the inset the In- 
dependence of the stretching parameter (3k is also re- 
ported. We remind that our model reproduces only the 
structural (a) relaxation processes usually explained in 
term of intra-basins transitions, at variance with other 
fast relaxation processes taking place inside the basins. 
This explain the presence of only one step relaxation in 
C(t). We also observe that: i) the relaxation dynamics 
is well represented by a stretched exponential decay; and 
ii) the stretching parameter (3k is strongly temperature 
dependent, implying a violation of the time-temperature 
superposition principle |9) . Moreover (3k decreases from 
1 at high T (Debye relaxation) down to w 0.35 at low 



T, a value consistent with experimental findings |1C| ] and 
theoretical prediction |ll[] in fragile glass-formers. 

In fig. 3a we show the shear viscosity and the relax- 
ation time r versus inverse temperature. They are almost 
proportional to each other, strongly increasing in a small 
temperature range (150 — 20 K). However, their tem- 
perature behaviour is well represented by an Arrhenius 
law, and does not show the dramatic increase expected 
for fragile glass-formers. Whether this unexpected be- 
haviour i) has to be ascribed to a failure of our model, 
or ii) is a genuine behaviour of LJ liquids at constant 
density, is still unknown. At those temperature where 
the direct MD calculation of n is affordable (o in fig. 3a) 
we found a good agreement between MD's and model's 
results, supporting the hypothesis ii). 

In fig. 3b we report the inverse diffusion coefficient, 
D~ x , versus the ratio rj/T. The Stokes-Einstein (SE) 
relation would predict direct proportionality, i. e. D oc 
T jr\. The full line (slope 1) indicates that, at high T, 
the SE relation asimptotically holds. Upon decreasing 
T, the slope £ decrease towards £ rs 0.28, indicating a 
breakdown of the SE relation, as observed in different 
experiments Jl2| . In particular, the crossover between 
the two regimes occours in the same temperature region 
where (3k deviates from 1. It is tempting to note that 
the crossover position, at -q/T « 0.1 Poise/K, and the 
fractional exponent at low T, £ ss 0.28, are in fairly good 
agreement with the experimental results in the fragile 
glass former o-terphenyl |K^jL4| ]. 

In conclusion, we presented a simplified model, based 
on a vibrational local dynamics and on collective jumps 
among minima, that well describes the structural relax- 
ation features of supercooled liquids and glasses. Explot- 
ing this model we are able to investigate the long time 
(low temperatures) dynamics. We recover, in the simple 
LJ system, some important features of real glass-former, 
in particular: a) stretching of the relaxational dynamcs; 
b ) failure of the time-temperature superposition principle 
(T-dependence of (3k)', and c) breakdown of the Stokes- 
Einstein relation. 

We wish to thank B. Coluzzi, G. Monaco and 
F. Sciortino for useful discussions, and D. Leporini who 
called our attention to the fractional SE relation. 
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Fig. 1 - Potential energy of the LJ system as determined from MD 
heating the glass (o) and cooling the liquid (•), and from the 
present model using all the minima (dotted line) or only the 
glassy minima (full line). 

Fig. 2 - Normalized autocorrelation functions of the off-diagonal 
elements of the stress tensor at the indicated temperature 
as determined from the model (open symbols). The lines 
represent the best fits to the data with a stretched exponen- 
tial time decay. The inset show the T dependence of the 
stretching parameter 

Fig. 3 - a) Shear viscosity and relaxation time versus inverse 
temperature, b) Inverse diffusion coefficient, -D , versus 
the ratio rj/T. According to the Stokes-Einstein relation, 
a linear proportionality is expected (dashed line). The full 
line, with slope 0.28, is the best fit to the low temperature 
data. 
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Fig. 3b - Connected network of minima.... - Angelani et al. 



